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Abstract 

Conceiving a molecnle as composed of smaller molecular fragments, or subunits, is one of the pillars 
of the chemical and physical sciences, and leads to productive methods in quantum chemistry. Using 
a fragmentation scheme, efficient algorithms can be proposed to address problems in the description of 
chemical bond formation and breaking. We present a formally exact time-dependent density-functional 
theory for the electronic dynamics of molecular fragments with variable number of electrons. This new 
formalism is an extension of previous work [Phys. Rev. Lett. Ill, 023001 (2013)]. We also introduce a 
stable density-inversion method that is applicable to time-dependent and ground-state density-functional 
theory, and their extensions, including those discussed in this work. 


1 Introduction 

Simple and productive methods to investigate dy¬ 
namical features of solids and molecules are of¬ 
fered by Time-dependent Density-functional Theory 
(TDDFT) [1]. TDDFT embodies many concepts and 
formal exact results, but its core is the 1-1 corre¬ 
spondence [2] between time-dependent (TD) external 
potentials and TD electronic densities, provided the 
initial state of the system is given. Every observable 
of the system is expressed as a TD density-functional. 
The calculation of the density is carried out by solv¬ 
ing the TD Kohn-Sham (KS) equations [3], which are 
single-particle Schrodinger equations that require an 
approximation to the exchange-correlation (XC) po¬ 
tential, a density-functional. 

The adiabatic local density approximation (ALDA) 
[4] to the TD XC potential is the simplest, use¬ 
ful approximation to study the dynamics of atoms 
and solids. However, when applied to molecules, 
ALDA often yields unphysical results, for example, 
atoms with fractional charges at dissociation, under¬ 
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estimated charge-transfer excitation energies, missing 
double excitations, among others. Two decades of re¬ 
search have shown that it is challenging to enhance 
the performance of ALDA while preserving computa¬ 
tional simplicity and elegance. 

The TD KS equations describe all of the electrons 
as part of a single entity, imposing a limit on the 
number of atoms that can be simulated in a reason¬ 
able amount of time. This limit can be increased 
by dividing a molecule into fragments and perform¬ 
ing inexpensive calculations on each individual frag¬ 
ment. Several approximated methods to investigate 
the electron dynamics of molecules as composed of 
fragments are available [5-8]. These methods typi¬ 
cally assign a set of TD single-particle Schrodinger 
equations (not necessarily TD KS equations) to ev¬ 
ery fragment in the molecule. The fragment elec¬ 
trons are subject to the usual interactions as if they 
were in the presence of the fragment nuclei only; and, 
an extra potential, usually fragment-dependent, ac¬ 
counting for the interaction between the fragments, 
is added. Successful applications to the calculation 
of solvachromatic shifts [9, 10] and excitation energy 
of monomers [6] have been reported. 
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A rigorous extension of TDDFT for molecules 
made of chemical fragments was presented in Ref. 

[11] , In this extension, a molecule is divided into 
fragments, each a set of atoms. Every fragment is as¬ 
signed an initial state, and a Hamiltonian including 
a global, auxiliary potential, termed partition poten¬ 
tial, which enforces that the total electronic density 
is the true TD electronic density of the molecule. We 
proved that the partition potential is uniquely deter¬ 
mined by the TD electronic density of the system, and 
that it can thus be expressed (and approximated) as a 
density-functional. The linear response and an exten¬ 
sion to consider electromagnetic fields are presented 
in Ref. [12]. The Hamiltonians used in Refs. [11] and 

[12] , and the aforementioned approximated methods, 
are particle-conserving, i.e., the average number of 
electrons in a fragment is time-independent; this re¬ 
striction is eliminated in this work. 

In this paper we first introduce a density-inversion 
method that improves that of Ref. [11] and allows 
us to estimate the partition potential. Using a model 
for the errors, we discuss, in close connection with 
standard TDDFT, the origin of uncertainties of the 
partition potential in regions far from the molecule. 
We extend our fragment-based TDDFT [11] to al¬ 
low for variable numbers of electrons in each frag¬ 
ment, while preserving the uniqueness of observables 
as density-functionals. The formalism introduced in 
this paper serves as a theoretical foundation for the 
development of methods accounting for electronic ex¬ 
citations and electron-transfer processes, without sac¬ 
rificing explicit use of the electronic density and com¬ 
putational efficiency. 

2 Fragment-based TDDFT 

2.1 Formulation 

An electron in a fragment, labeled a, is subject to a 
1-body external potential, denoted as Va- For exam¬ 
ple, Uo,(r) = “ R-d; la is a set of labels 

for the nuclei in fragment a. We assign to each frag¬ 
ment in the molecule a Hamiltonian that includes an 
auxiliary potential, the partition potential Vp'. 

Hc[vp]{t) = H°+ J d^r n{r)vp{rt) , (1) 

where = T -\-W f d^r fi(r)uQ,(r), T and W are 
the kinetic, and coulombic repulsion energy opera¬ 
tors, respectively, and n(r) is the density operator. 
This Hamiltonian does not include external driving 


forces other than those due to the nuclei of fragment 
a. TD displacement of the positions of the nuclei can 
be described by introducing a time-dependent Hamil¬ 
tonian where Va is replaced by the corresponding TD 
fragment-potential, 

The state of a fragment is described by the evolu¬ 
tion of the ket |V'akp](^)) Fock space, which sat¬ 
isfies the TD Schrodinger equation (atomic units are 
used throughout): 

idt\-lpa[Vp]{t)) = Ha[Vp]{t)\'lpa[Vp]it)) , (2) 

where 

li’ait)) ■ ( 3 ) 

M 

{V'a. m } are kets corresponding to states with integer 
number of particles and {h'a,M} are the weight am¬ 
plitudes of those states. Kets with different number 
of electrons are orthogonal, (V'a.MjV'a.M') = 0 , M ^ 
M'. The total density is given as 

n{rt) ='^na{rt) , (4) 

a 

which defines Up when u„(rt) = {ipa{t)\n{r)\'tljc,{t)) , 
as proven in Ref. [11]: Given a set {V'a, 0 ) "ya}) two 
potentials Up and v'^ that differ by more than a TD 
constant cannot give rise to the same density. A 
corollary of this theorem is that there is a TD density- 
functional that, when evaluated at a given TD elec¬ 
tronic density, gives the corresponding TD partition 
potential 

The partition potential represents the TD elec¬ 
tronic density of the supermolecule, and is decom¬ 
posed as follows [12]: 

Vp{rt) = VG{rt) -I- Vd{rt) . (5) 

VQ is a “gluing potential”, accounting for the corre¬ 
lation between the fragments, and Vd is the driving 
potential the molecule is subject to (e.g. laser field). 
The gluing potential yields the shape of the potential 
such that the TD electronic density is recovered. The 
gluing potential satisfies [12]: 

yV •u(rt)VuG(rt) = {'tlj(t)\[H° ,V ■ j{r)]\tp{t)) 

-Xl^^“k)l[Ra.V-j(r)]|V'a(0) • ^ 

oc 

^This theorem and corollary were recently used in Ref. [13] 
to propose an inversion method in the context of embedding 
potential-functional theory. 
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The terms on the right-hand-side of Eq.(6) are TD 
density-functionals. Approximating these terms and 
solving the resulting differential equation yields an 
estimate of the gluing potential. Another route to 
approximating vq is by assuming that the system 
evolves adiabatically from its ground states, driven by 
a very slowly-varying field. In such case, the poten¬ 
tial VQ is obtained from the adiabatic approximation 
in ground-state Partition DPT [14, 15]: 

= , (7) 

where is the external perturbation the in¬ 

teracting electrons are subject to in their ground- 
state in order to yield the density n(rt) (the unique¬ 
ness of follows from the Hohenberg-Kohn theo¬ 
rem). The partition potential, t!^'^[n(t)], is the La¬ 
grange multiplier required to solve the minimization: 


min 




( 8 ) 


under the constraint of Eq.(4). The Lagrange mul¬ 
tiplier for this problem is unique, up to an arbitrary 
constant [16]. 

The TD partition KS equations are: 

-b vii^c[na]{r,t) 

^2 ( 9 ) 

-b Ua(r) -b VG[n]{r,t) + Vd{r, . 

The density can be obtained by means of: n(r, t) = 
, where {fia} are the (time- 
independent) occupation numbers, chosen from a 
proper ensemble [11]. 


3 Classical Interpretation of 
the Partition Potential 

Consider a system composed of a single massive par¬ 
ticle, and a bath made of particles much smaller 
than the massive one. All the particles in the par- 
ticle-bbath system interact via a potential, Umt, of 
the form (ri —^j)- The evolution of the sub¬ 

system particle, labeled S, is dictated by Eq. (2). 
There is no partitioning of the external potential be¬ 
cause the particle and the bath are subject, implic¬ 
itly, to a macroscopic, external, confining potential. 
The Hamiltonian of the subsystem particle is thus 
Hs = T + J d^r h{r)vp{rt), and the Hamiltonian of 
the bath is Ab = T-b IE-b/ d^r h(r)up(rt). The aver¬ 
age position of the particle is fs(t) = / d^r r ns{r,t), 
where ns(r,t) = \'tp{r,t)\‘^. 


By the Ehrenfest theorem and correspondence 
principle we have 

=Fp^sit) , ( 10 ) 

where Fp_s(t) = —f d^r ns(r,t)Vvp(r,t), and ms 
is the mass of the particle. Comparison with 
the equation of motion of the real system indi¬ 
cates that, in the classical limit, (Vvp)(fs(t)) = 
(VrsUint)(rs(i),rB(t)), where tb is the average co¬ 
ordinate of the bath. The partition force at the posi¬ 
tion of the particle is simply the force exerted on the 
particle by the bath. 

As the mass of the subsystem particle is increased, 
the density tends to a Dirac distribution. It follows 
from the above result and Eq. (6) that the shape 
of the partition potential for any point but that of 
the particles is indefinite. However, for given ini¬ 
tial momenta and coordinates of the particles and 
bath (this condition replaces the requirement that 
the initial wavefunction be specified), the trajectory 
of the momenta of the total system is in a one-to- 
one correspondence with the trajectory of partition 
forces exerted on each particle. Furthermore, if the 
assumptions of Langevin dynamics are applicable, 
the partition force on the massive particle can be 
interpreted as Fp^s(0 = ~7Vs(i) + Fran(i)) where 
vs(t) = dfs/dt, 7 is the friction coefficient, and Fran 
is the random force. 

The gradient of Vp can only be known at the po¬ 
sition of the particles, nowhere else. In the origi¬ 
nal proof by Runge and Gross [2] , and its extension 
for partition potentials [11], it was found that poten¬ 
tials whose gradients grow rapidly in regions distant 
from the molecule are not covered. This result, as 
illustrated in Section 4, is related to the uncertainty 
in the estimation of the partition potential near the 
boundaries of the simulation box. 


4 Numerial TD Potentials 

TDDFT, in which our formulation is built upon, con¬ 
cerns itself with the simplification of the problem: 

(idt - II^lv](t))ji/;(t)) = 0, IV'(O)) = IV'o) , (H) 

where 

H^[v]{t) = T + \W + J d^r n{r)v{rt) . (12) 

Runge and Gross [2] showed that if v is Taylor- 
expandable around t = 0 and does not have physi¬ 
cal anomalies in the boundaries, then v determines n 


3 



uniquely up to a TD constant in the potential; this 
theorem (recently used to formulate quantum con¬ 
trol in TDDFT [17]) can be extended to include non- 
analytic potentials [18]. Let us denote the RG map 
as thus, n{t) = A^^[ri](t). The operator W can 
stand for different types of electron-electron interac¬ 
tions, such as screened coulombic repulsion. If A = 0, 
then the electrons are non-interacting. 

Suppose a well behaved density, and an ini¬ 
tial state, tpOi are known. If vi and vq exist, where 
VA(i) = [n‘'®^](t), then the Hartree-exchange- 

correlation potential for the system, by definition, is 
vhxc = Vo — Vi- Using the exact map would re¬ 
quire solving the TD Schrodinger equation, which is 
what one wants to avoid in practical calculations. 

For the development of functionals, exploration of 
the map A^^ is fruitful; this map could be investi¬ 
gated by solving the problem — A^^[u](t) = 0, 

which is a root-hnding problem. The first order re¬ 
sponse of the density for some perturbation Sv is 
Sn{rt) = J d^rdF rT')i5u(r't'). The response 

function x~^ should decay in the asymptotic regions 
because large perturbations of v in those regions have 
a small response in n. This relation between densities 
and potentials introduces instabilities in root-finding 
algorithms aimed at reproducing the density corre¬ 
sponding to a given external potential. In the ground- 
state case, the instabilities could be eliminated by 
enforcing satisfaction of eigenvalue constraints. For 
three dimensional applications, in general, capturing 
the asymptotic regions is difficult when Gaussian ba¬ 
sis sets are used because they do not have the correct 
asymptotic behavior. 

Instead of looking for the exact root, one can solve 
a minimization problem: 

min / l|n"®^(s) - A^ [u](s)ll2 ds , (13) 

where jj/jj is a suitable spatial norm, and V is a 
space of physical potentials. The quantities n''®^(t) 
and {'ip[v]{t)\n{r)\'ip[v]{t)) need to be approximated. 
Let us write n"'®^ (t) — A^^ [ti] (t) = h"'®^ (t) — A^^ [?;] (t) + 
e[n''®^z;]. fi''®^(f) is the approximation to n’'®^(t) and 
A^j, b] is the approximation to A^^. If v* is the ex¬ 
act potential representing n’’®^, then the problem be¬ 
comes h“'®^ = A^^ b*] -I- e. Because we cannot use 
exact methods to determine n“'®^ and A^^, we assume 
that e is a random function of the space-time coordi¬ 
nates. Moreover, one would expect that h“'®^ and A^^ 
have smooth timespace gradients, and that e displays 


autocorrelation because the spacing between points 
is arbitrarily small. 

4.1 Estimation of the Partition Poten¬ 
tial 

Let Vp be a space of TD partition potentials, and V 
a space of TD densities and define the map: 

Asg-.Vp^V, (14) 

where Sq = {ipa, 0 jVa}- For a given TD partition 
potential, the density is obtained by evaluation of the 
above map at the given partition potential; in other 
words, n{t) = AsQbp](f). This map depends on the 
history of the partition potential, i.e., it has memory 
dependence [11]. 

Let V* be the true partition potential. We assume 
that, due to numerical errors, the estimation to the 
reference density h“'®^ is of the form tt,’’®^ = Asg -f e, 
where e is a random function. To estimate the parti¬ 
tion potential corresponding to h“^®^ we minimize: 

lbbp]llM = lb’''^'-Asobp]ll^ (15) 

where is the measure. 

Since e is a function, its probability density func¬ 
tion (PDF) is a functional. The PDF depends on 
parameters, denoted collectively as 0, and the PDF 
itself as D([e]]0). The probability that e is observed 
in a set If is given by the path integral: 

P{e e U\Q) = [ dmL[e] D([e]10) , (16) 

Ju 

where the measure over the space of errors is tol. 
The traditional methods of non-linear regression can 
be applied to estimate the best parameters of the dis¬ 
tribution, 0*, for a given set of observations. Then a 
Taylor expansion in terms of the parameters can be 
used to generate their PDF, which can then be used 
to estimate the error in the parameters. In this case, 
the parameters are the variance and the partition po¬ 
tential. 

In the next section, we will expand the partition 
potential in a spline basis set. The parameters are 
the values of the partition potential at the knots, and 
they are correlated: A perturbation of the partition 
potential at one knot affects the response of the den¬ 
sity in other knots. Hence, we must employ a model 
of correlated errors. Finding the correct model is a 
demanding task beyond the scope of this work. For 
this reason, we choose a biased model based on the 
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following observations: i) A measure of the error of 
the form / d^rdt t) — A 5 (j[up])^ suffers of au¬ 

tocorrelation. ii) Far from the molecule, the parti¬ 
tion potential has little influence on the density, iii) 
Estimating the density is not sufficient; its spatio- 
temporal gradient is an important quantity. An error 
measure accounting for these observations is: 
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l|ebp]||^ = / + (9te(r,f))2} . 

(17) 

Based on ii), we choose a measure of the form 
d^(r,t) = d^rdt ^t)(5(r — r^). Where {r^} 
are points selected in such a way that |Vep -I- (9*6)^ 
resembles a y^-distribution. To apply this measure 
of error in the next section, we transform the above 
measure into a vector norm. Then, the resultant dis¬ 
tribution is expanded in terms of the partition poten¬ 
tial evaluated at the knots, and asymptotic analysis 
is applied [19], leading to the random variables re¬ 
quired to reproduce the density within a small error 
tolerance. 




4.2 Id Electron in a Double-well Po¬ 
tential 

Let us revisit the example considered in Ref. [11]: A 
one dimensional “electron” in a double well potential. 
The TD partition equations are: 

= (^ - -dl + Va{x) + Vp{x,t)j4>a{x,t) , 

(18) 

where a = L, R, standing for left and right wells. 
The potentials are Va (x) = vq/{x — Xa)^ + a; the 
parameters are: vq = —1, xr — a:L = 4, and a = 1. 
The density is obtained by averaging over the orbital 
densities of both wells: 

n{x,t) = ^\<j)L{x,t)\'^ + ^\(j)Rix,t)\'^ . (19) 

Suppose that the supermolecule evolves from the 
ground state driven by a monochromatic laser. The 
evolution of the system is thus dictated by the solu¬ 
tion of: 

idttpix, t)= + v{x) + Vd{x, t)'^tp{x, t) , (20) 

where Vd{xt) = Ex sinuit^ and the external potential 
is u = ul -b Ur. The density obtained from the above 
evolution equation is n“'®^(xt) = |'!/'(xt)l^, which is the 
target density we wish to represent. 


Figure 1: Snapshots of the partition potential. In 
a), solid line: Initial partition potential, dashed line: 
Left fragment external potential, dashed-dotted line: 
Right fragment external potential. In b), d), and 
f), solid lines: Left electronic density, dashed lines: 
Right electronic density. In d) and f) the dashed- 
dotted line is the total density. 


The laser parameters are uj = 0.3, Eq = 0.05. We 
propagate the states of the system using the Crank- 
Nicholson method; time step is 0.1, box length is 
20, spatial step is 0.2, and total propagation time 
is 10 units. The partition potential is represented 
in a smoothed, cubic, spline basis set with 22 knots 
equally spaced in the box. The initial partition po¬ 
tential is estimated by minimizing the error using se¬ 
quential quadratic programming (other useful meth¬ 
ods employ the density error directly [20]). To obtain 
the initial densities, the error functional of Eq. (17) 
(with the time-dependency dropped) is minimized. 
First the problem {-l/2dl + Va+v°)cj)n,a = en,a4‘n,a 
is solved (with the finite-difference method) for both 
wells with some estimation of v^. Then, the den¬ 
sity is compared with that of the system of reference 
in order to propose the next estimation in the itera¬ 
tive procedure of sequential quadratic programming; 
the constraints are conservation of charge and chem¬ 
ical potential (HOMO) equalization [12]. Figure l.a. 
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shows the initial partition potential and external po¬ 
tentials of each well. The initial fragment densities 
that add up to the ground-state density of the super¬ 
molecule are displayed in Figure l.b. 

The estimation of the evolution of Vp is determined 
by using the step-by-step scheme proposed in Ref. 
[11], and the norm of Eq. (17): For example, for 
At = 0.1, the error norm can be approximated as 

At At)[{dxe{xi, At))^ -I- {dte{xi, At))^]} . 

( 21 ) 

The numerical value of the above quantity depends on 
the value of Vp at the spline knots and at t = Atj2. 
Hence, this quantity is varied until the above func¬ 
tion is minimized and the total density of the frag¬ 
ments match that of the true system. The procedure 
is repeated similarly for the rest of the propagation. 
Figure l.c. shows the partition potential at t = 1.0; 
it is localized in the intermediate region between the 
fragments. The fragment electronic densities (Fig¬ 
ure l.d) are also well localized at t = 1.0. Because 
in absence of the partition potential the fragment- 
densities would just be localized around their wells, 
the partition potential must be such that it induces 
the transfer of charge from the right fragment into 
the left fragment (Figure l.e). However, as we note 
in Figure l.f, the charge transfer in this case is repre¬ 
sented by the spreading of the right electronic density 
into the left one, not by a change in the fragment pop¬ 
ulations. Two observations: i) If one were to assign 
a grid that is fine enough around the center of the 
wells and then coarser as one moves away from the 
wells, then to describe the density spreading, the grid 
would need to be time-dependent, ii) The partition 
potential must induce the charge transfer and act like 
a “spoon”. 



Figure 2: Error-bar plot of the partition potential at 
t = 6.2. Solid line: Interpolated optimal potential, 
circles: Interpolation knots. 


The result of the error estimation in the partition 
potential at t = 6.2 is shown in Figure 2. As ex¬ 
pected, in the boundary regions, the error is quite 
significant, and the error bars extend well beyond 
the plotting range. This implies that the shape of 
the potential in these regions is not reliable. All 
space-time points obeying causality are coupled. For 
example, variation of a single knot in the boundary 
affects its neighbors, introducing large gradient fluc¬ 
tuations. Thus, the error can spread to regions where 
the density is non-negligible. This can cause instabil¬ 
ities in the minimization procedure if a norm such 
as / d^rdf (^^^(r,^) — n(r,t))^ is employed. For this 
reason we recommend the use of derivatives of the 
density to measure the error. 

5 Variable Occupation Num¬ 
bers 

Let us assign variable electron-occupation numbers 
to the fragments. First, divide the total propagation 
time into blocks:[0, r) U [r, 2t) U ... U [{m — 1)t, mr), 
where mr is the final time of the propagation, and 
let 

= (22) 

be a set of instantaneous kets, in Fock space, for frag¬ 
ment a. At a single time t = kr, the following min¬ 
imization is performed to obtain the set of kets de¬ 
scribing the density of the fragmented molecule: 

{ie^)}^lT =argmin{^(V^S|iJ°|V^S) s.t. 

oc 

=n{r,kT), Vrj , 

O' 

(23) 

The occupation numbers of fragment a are for¬ 
mally expressed as \va,M{kT)\'^ = |(?a,Ml?a)P • Here, 
is an optimal ket (obtained from solving the 
right hand side of Eq. (23)) for fragment a with M 
electrons. These numbers are density-functionals. 

The evolution operator of fragment a is: 
Ua[vp\{ti,t2) = Texp(-ids I7 q,[up](s)) . Intro¬ 
duce the displaced set of kets: 

A„ = {t/„(r,0)|C°), 

17„(2t,t)|^^), ... UaijriT, (m - I)t)|C"^)} • 

(24) 

Now let us define the following dyadic product: 
{X^Xl){k) = The symbol A„At is the 
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set of dyadic products where the fc-th component is 
the dyadic product between the ket at the beginning 
of the fc-th block and the displaced ket from the k — 1- 
th block. Now, let Ba be the TD operator: 

k=l 

where Wr is the Dirac-Comb kernel. Addition of the 
operator Ba to the Hamiltonian Ha (t) yields the non- 
Hermitian operator: 

-ffc.abp](f) = Ha[Vp]{t) + iBa[Vp]{t) . (26) 

The evolution of the system is now determined by 
IV'abp]), which obeys 

i^t\^pa[Vp]it)) = Hc,a['>^p]{t)\tlJa[Vp]it)) . (27) 

The total density is n(rt) = 

and the number of particles in fragment a is Na (t) = 
{'tpait)\N\ipa{t))- In general, any observable, 0{t), is 
expressed as a functional of the partition potential, 
{i^a[Vp]{t)\d{t)\lpa[Vp]). 

Given the partition potential and occupation num¬ 
bers as density-functionals, the scheme to determine 
the evolution of the molecule is the following: First, 
propagate the kets Htpa)} in the interval [0,r) with 
fixed populations on each fragment. Then, at t = r, 
obtain new occupation numbers according to Eq. 
(23) as well as new states to propagate; continue 
the propagation in the block [r, 2r). The proce¬ 
dure continues similarly for the rest of the propa¬ 
gation. The density of the system is then obtained 
as n{r,t) = X]a(V’a(O|n(r)|'0a(O)- The theorem dis¬ 
cussed in section 2 also applies in this case. Therefore, 
the partition potential for this scheme is uniquely de¬ 
termined by the TD electronic density, up to an ar¬ 
bitrary constant. 

The partition potential is discontinuous at the re¬ 
laxation nodes (points where t is an integer multi¬ 
ple of r). Discontinuities in time can be eliminated 
by using an integral transformation that smooths the 
observable at the relaxation nodes. In practice, how¬ 
ever, it is convenient to propagate the occupation 
numbers and gluing potential assuming that they are 
continuously differentiable functions of time. It can 
be shown, assuming that the dynamics of the occupa¬ 
tion numbers is much slower than that of the partition 
potential, that the I-I map between the former and 
the density still holds. This follows from the scheme 
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Figure 3: Evolution of the fragments with TD elec¬ 
tron populations. In a) the solid line is the result from 
the inversion, the dashed line is the result from the 
two-state approximation, and the dashed-dotted line 
is obtained by omitting the gluing potential: vq = 0. 
In c) and e), solid line: til, dashed line: ur, dashed- 
dotted: n. 


we have shown here because the electronic popula¬ 
tions are fixed in the first block, allowing us to apply 
the Runge-Gross theorem in such block. 

A density-functional approximation to the occupa¬ 
tion numbers is needed to apply the theory. The 
dynamics of the occupation numbers can be inves¬ 
tigated using master equations, where the rate co¬ 
efficients are determined by Eermi’s golden rule, or 
transition elements that couple the fragments. Here, 
we illustrate a simple approach: A trial wave function 
to investigate the evolution of the occupation num¬ 
bers is \r]{t)) = aL(t)|<PL) + aR(f)|(/3R), where \ipa) 
is the ground-state of the electron described only by 
Ha (This hamiltonian in coordinate representation is 
— 1/23^ -I- '(;Q,(a;)). The dynamics of electron trans¬ 
fer is governed by a two-component wave-function 
a = (ol, aR)"*". We assume that the Hamiltonian cou¬ 
pling that relates the two fragments is of the form: 

Hit) = Hf + f dx {vg{x,0) + Vd{x,t))n{x) (28) 
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where is the uncoupled Hamiltonian; 

Hal^Pp) = 0 if a ^ /3. For the sake of the illustration, 
the gluing field is frozen, serving as a “bridge” for the 
charge to be transferred from one well into the other. 

From the evolution equation: idt\r]{t)) ='H{t)\ri{t)) 
we infer that the state vector, a, satisfies: 

idtSL{t) = S“^(eo + A(t))a(t) (29) 

where Sajs = / dec (p*(x)(p/3(ce), co = diag(eo,eo), and 

Aa/s(t) = / dx lfi*(x)(vG(x,0)+Vd(x,t))lfi/3(x) . 

(30) 

The occupation numbers are obtained from the “den¬ 
sity” of a: Na(t) = laal^(t) + Re(al(t)aR(t)SLR). 
The last term arises from the overlap of the functions 
(fiL and (pR, guaranteeing that A^l + A^r = 1- 

The example of the previous section, summarized 
in Fig. 1., illustrates how the partition potential can 
be estimated even under a strong laser field. We now 
return to the example of section 4.2 and allow for 
variable occupations. The parameters for the propa¬ 
gation now are t = 2, At = 1, uj = 0.3, i?o = 0.02. 
Comparison of the exact time-dependency of the av¬ 
erage number of electrons of the left fragment with 
the approximation described above is shown in Fig¬ 
ure 3.a. The initial gluing potential used here is that 
shown in Fig. l.a. The two-state approximation 
works well at short times, and displays deviations af¬ 
ter t = 20. Capturing the results of the two-state ap¬ 
proximation would be quite challenging by fixing the 
occupation numbers and finding the corresponding 
partition potential. Improvements over the two-state 
approximation can proceed by either refining the glu¬ 
ing potential (going beyond the frozen approxima¬ 
tion) or increasing the number of states considered 
to couple the fragments. The first alternative has the 
advantage that the equations can be solved very fast. 
Nonetheless, for functional development, the gluing 
potential is also a determinant factor for the evolu¬ 
tion of the shape of the electronic fragment-density 

(/ Ua/Na). 

Figure 3.b shows a snapshot of the “exact” par¬ 
tition potential at t = 40. In contrast with the re¬ 
sults of section 3.2, the partition potential is now well 
localized (Figures 3.d and 3.f). This suggests that 
the standard methods of ground-state PDFT can be 
used to estimate the partition potential through the 
use of the adiabatic approximation. The fragment 
densities also remain localized (Figure 3.c, 3.e, and 
3.g). Qualitatively, the partition potential accounts 


for the shape of the electronic densities of the frag¬ 
ments, while the occupation numbers are responsible 
for their height. 

6 Concluding Remarks 

We formulated a TDDFT for treating molecules as 
composed of smaller composite units. The successful 
application of this formulation requires approxima¬ 
tions to the partition potential and the occupation 
numbers. This can be accomplished by a proper ap¬ 
proximation to the Hamiltonians {Hc^a(t)}, or the 
auxiliary evolution equations of the electron popula¬ 
tions in the fragments. The error analysis was also 
presented. It leads to a simple form of estimating 
the errors in the potentials. The partition potential, 
obtained by numerical inversion, can be uncertain in 
spatio-temporal regions where the density is small. 
However, as time increases, the error can propagate 
from the boundary areas into regions were the density 
is high. 
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